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We numerically investigate the competition between phase separation and dynamical arrest in a 
colloidal system interacting via a short ranged attractive potential. Equilibrium fluid configurations 
are quenched at two different temperatures below the critical temperature and followed during 
their time evolution. At the lowest studied T, the phase separation process is interrupted by the 
formation of an attractive glass in the dense phase. At the higher T, no arrest is observed and 
the phase separation process proceeds endless in the simulated time window. The final structure 
of the glass retains memory of the interrupted phase separation process in the form of a frozen 
spinodal decomposition peak, whose location and amplitude is controlled by the average packing 
fraction. We also discuss the time evolution of the non ergodicity parameter, providing evidence of 
a progressively decreasing localization length on increasing the packing fraction. Finally, we confirm 
that the reported results are independent on the microscopic dynamics. 



PACS numbers: 61.20.Lc, 64.70.Pf, 47.50.+d 



I. INTRODUCTION 

Advances in colloidal science make it possible to real- 
ize systems of interacting particles, with a tunable inter 
particle potential Interactions can be controlled in 
the range and in the strength, expanding considerably 
the possibilities found in atomic and molecular systems. 
The size of the colloidal particles allows experimentalist 
to study these systems with a wide range of optical tech- 
niques, such as Dynamic Light Scattering and Confocal 
Microscopy. Among the colloidal systems which have no 
atomic counterpart are the so-called short ranged attrac- 
tive colloidal (SRAC) systems, in which the inter-particle 
potential has a range significantly smaller than the col- 
loidal size. Experimentally these system can be realized 
by adding to a solution of colloidal particles a depletant 
agent, generally a polymer, When the radius of gyration 
is smaller than the colloidal size, the large particles ex- 
perience an effective attraction whose range is related to 
the size of the polymer and the intensity to their concen- 
tration 0. 

SRAC systems posses an extremely rich thermodynamic 
and dynamic behavior (for recent reviews see for example 
Refs. 0, 0, 0, 0)- The short range attraction affects pro- 
foundly the structure of the phase diagram. When the 
range of the attraction is comparable to the particle size, 
the phase diagram presents a typical structure of a van 
der Waals fluid, i.e. a liquid-liquid critical point at low 
density and a solid-fluid phase transition at high density. 
However, when the range is much shorter than diameter 
of the colloidal particles, the liquid-liquid phase separa- 
tion becomes metastable and is buried inside the fiuid- 
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solid coexistence curve. This phenomenompredicted the- 
oretically within perturbation theory LP, H, B ■ ha s been 
confirmed by simulation jlOl ] and experiments [ill Il2j . It 
is interesting to note that critical fluctuations related to 
metastable critical point favor nucleation of the crystal 
phase |l3j. 

Slow dynamic properties in SRAC systems also show fea- 
tures not observed in atomic systems. When short ranged 
attractions are present, the fluid-glass line (as well as its 
precursors, the iso-diffusivity lines 0,0]) shows a reen- 
trant behavior in the interaction strength-packing frac- 
tion plane. Two different competing mechanism for dy- 
namic arrest are present, generating respectively a repul- 
sive dominated and an attractive dominated glass. Dy- 
namics become so complex that the usual stretched ex- 
ponential decay characterizing the usual slowing down 
close to a a glass transition, cross over toward a loga- 
rithmic decay [jj, [H E M, H E| • These novel phe- 
nomena, first proposed on the basis of theoretical cal- 
culations based 0, on the Mode Coupling Theory 
developed by Goetze and collaborators J23| have been 
subsequently confirmed by simulations 0, 0, [23|, 
and experiments pi IM I27ll2c| . 

An additional phenomenon which is not usually encoun- 
tered in atomic systems takes place at low packing frac- 
tions, where SRAC systems are known to form gel, i.e. 
arrested space spanning structure 0, [3(1 El E3- 
The nature of the gel transition in short-range attractive 
colloidal systems has received significant attention in re- 
cent years (for a recent review see for example Ref. HH). 
Several routes to the gel state have been proposed and 
critically examined, with a special emphasis on the analo- 
gies and differences between glass and gel formation. It 
has been shown that, for some class of potentials, struc- 
tural arrest at low packing fractions can be interpreted as 
a glass phenomenon of self-assembled clusters [3a . l3rjl . l37j . 
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For the case of SRAC, evidence is building that the gel 
state results from an arrested phase separation when the 
phase separation dynamics generate regions of local den- 
sity sufficiently large to undergo an attractive glass tran- 
sition [H HI E3, El El El El . 

In this manuscript we report results of Newtonian and 
Brownian molecular dynamics simulation of a binary 
mixture of particles interacting via a short range attrac- 
tive square well. For this system, the location of the 
phase separation region and of the repulsive and attrac- 
tive lines have been previously studied. The attractive 
glass line intersect the colloidal-rich colloidal-poor coex- 
istence curve on the colloidal rich side at an intersection 
temperature T x , a few per cent smaller than the critical 
temperature T c . We quench initial configurations at dif- 
ferent densities, equilibrated at a temperature where the 
system is homogeneous (T >> T c ) to two temperatures 
lower then the critical one (T < T c ) and study the coars- 
ening dynamics. For the case of the higher final tempera- 
ture, the coarsening dynamics never arrests and the sys- 
tems slowly approached to equilibrium phase separated 
state characteristic of a conserved order parameter phase 
separation process. For lower temperture quench, we find 
that the coarsening dynamics arrests within the simu- 
lated time, leaving the system in an out-of-equilibrium 
arrested structure. The morphology of the resulting gel 
is controlled by the phase separation process at large dis- 
tances (or small wave- vectors) . At length scales compa- 
rable to the particle size, the structure resembles the one 
characteristic of the attractive glass. 

II. DETAILS OF THE NUMERICAL 
SIMULATIONS 

We investigate a system that has been extensively 
studied earlier, a binary square well (SW) mixture |l5l 
HJ. The binary system is a 50%-50% mixture of N = 
2000 particles. The two species (labeled A and B) are 
characterized by a diameter ratio oaI^b = 1.2, which 
effectively suppress crystallization. Masses are chosen to 
be equal and unitary, i.e. m a — mb — 1. The SW inter- 
action is defined according to: 

!oo r < a a £ 
-U <J a ,p < r < Cr Qi(3 + A Qj/3 (I) 
r > <7 a .j3 + A a ^ 

where a a ^ — (a a + <rp)/2, a, (3 — A, B and A aj| g is 
the range of the attraction. We fix a a> p and the well- 
width e = . A ° ,/:l — = 0.005. The chosen e value is 
arbitrary, but representative of all SW potential with 
interaction range smaller than a few percent as far as 
thermodynamic and dynamics equilibrium properties are 
concerned El EH- For the chosen potential T c » 0.20. 
We choose fc# = 1 and set the depth of the potential 
Uq = 1. Hence T = 1 corresponds to a thermal en- 
ergy kgT equal to the attractive well depth. The diam- 
eter of the small specie is chosen as unity of length, i.e. 



erg = 1. Density is expressed in term of packing fraction 
<j> = (pa&a + Pb&b) ' ^/G, where p a = N a /L 3 , L being 
the box size and N a the number of particles of specie 
a. Time is measured in units of erg • (m/uo) 1 ^ 2 ■ ND 
has been coded via a standard event driven algorithm, 
commonly used for particles interacting with step-wise 
potentials [46|. Between collisions, particles move along 
straight lines with constant velocities. When the dis- 
tance between the particles becomes equal to the dis- 
tance where the potential has a discontinuity, the veloc- 
ities of the interacting particles instantaneously change. 
The algorithm calculates the shortest collision time in the 
system and propagate the trajectory from one collision 
to the next one. Calculations of the next collision time 
are optimized by dividing the system in small subsys- 
tems, so that collision times are computed only between 
particles in the neighboring subsystems. BD has been 
implemented via the position Langevin equation: 

**(*) = ir!f + (2) 

coding the algorithm developed by Strating |47|. In Eq.[3 
Ti(t) is the position of particle i, fi(t) is the total force 
acting on the particle, Do is the short-time (bare) diffu- 
sion coefficient, r$(t) a random thermal noise satisfying 

< rj(t)rj(Q) >= kBT5{t). In Strating's algorithm, a ran- 
dom velocity (extracted from a Gaussian distribution of 
variance y/ksT/m) is assigned to each particle and the 
system is propagated for a finite time-step 2 f^- , accord- 
ing to event-driven dynamics. We chose D such that 
short time motion is diffusive over distances smaller than 
the well width. 

Initial configurations were equilibrated at T = I for den- 
sities ranging from <\> = 0.01 to <f> = 0.50. For each den- 
sity, the system was quenched at two different final tem- 
peratures, T f = 0.05 << T x and T f ~ 0.15 ~ T x . The 
constant T evolution was then followed in time. In the 
case of ND, the characteristic time of the thermostat as 
been chosen so that velocities reach thermal equilibrium 
before the system starts to rearrange its structural de- 
grees of freedom. 

III. RESULTS AND DISCUSSION 

A. Phase diagram 

The coexistence curve in the limit of vanishing well 
width, i.e. the celebrated Baxter limit |48j, has been 
recently precisely calculated by Miller and Frenkel E3- 
These data provides an accurate estimate of the coexis- 
tence curve of all short range interacting potential, via 
an appropriate law of corresponding-states E3l- Indeed, 
Noro and Frenkel noticed that, when the attraction range 
is much smaller than a few percent of the particles diam- 
eter, the value of the first virial coefficient at the critical 
temperature was essentially independent on the details of 
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the potential, down to the Baxter limit. Following this 
prescription we reported in Fig. ^ the coexistence curve 
calculated by Miller and Frenkel in the r-<j) plane (r is 
the Baxter stickiness parameter) after transforming it to 
the T-<fi plane according to 

1-^= [l-{e (3ua -1)[(1- e)- 3 -!]] (3) 

i.e. imposing the equivalence of the virial coefficient of 
the Baxter and binary SW model. To independently sup- 
port the mapping, we calculate explicitly the location of 
the spinodal line, by bracketing it with the highest T 
at which no phase separation takes place and the lower 
T at which we observe a growing spinodal decomposi- 
tion peak. The calculated spinodal line, also reported in 
Fig. ^ is in very good agreement with Miller and Frenkel 
estimates. 

Fig. 2] also shows isodiffusivity lines for the studied SW 
model, i.e. lines along which the normalized diffusion 
coefficient is constant These lines show the typical 
non monotonic behavior characteristic of SRAC systems, 
in which a speed up of the dynamics takes place at the 
temperature at which the competition between the two 
different arrest mechanisms is balanced. The isodiffusiv- 
ity lines provide also a guide to the shape of the glass 
transition line, which has been numerically estimated, 
by extrapolation, as the zero isodiffusivity line (24|- As 
previously reported, the glass line meet the coexistence 
and spinodal line on the colloid-rich side, confirming that 
arrest in SRAC systems can take place at low packing 
fractions only as a result of an interrupted phase separa- 
tion mim. 

Fig. 2] also shows the location of the studied points in the 
phase diagram. 

B. Potential Energy - Number of Bonds 

One of the advantage of the SW models is represented 
by the fact that the potential energy per particle U/N 
is directly related to the number of bonds rib per par- 
ticle by rib = —2U/(Nuo). Indeed two particle can be 
unambiguously considered bonded if their relative dis- 
tance is within the attractive well distance. In the initial 
configurations, i.e. T = 1.0, the potential energy varies 
between values characteristic of few particle bonded, for 
the lowest density, to less then one bond per particle, for 
= 0.50. 

Fig. [2^ reports the time evolution of U/N for the ND 
case and Tf = 0.05. Following the quench at t = 0, the 
velocity degrees of freedom equilibrate and before t ~ 1.0 
the selected Tf is reached. Structural degrees of freedom, 
however, relax on a much longer time scale, dragged by 
the tendency of the system to phase separate. The evolu- 
tion of the energy shows three different processes. After 
a characteristic time, which is longer the smaller the ini- 
tial packing fraction, aggregation sets in and the energy 
decreases significantly. This aggregation process slows 



down significantly once a number of bonds of the order 
of 6 per particle is reached. In this late stage of the simu- 
lations run, two different behaviors can be distinguished. 
At very low density, i.e . <f> < 0.05, the energy continues 
to drop and does not seems to reach any stationary value 
in the simulation time window. For (f> > 0.05 , however, 
the time dependence of the energy abruptly stops and 
the system does not show sign of further evolution. This 
can be seen as first indication that an arrested structure 
has formed. 

Fig. [2t> reports data similar to the one presented in 
Fig.|2Ji but for BD. Despite the different microscopic dy- 
namic, the time dependence of the energy shows similar 
trends and the value of the energy at which the system 
stops for different densities is indeed very similar to the 
ND case (see Fig. |3J). Comparing in more details the 
two dynamics, one notice that the Brownian aggrega- 
tion dynamics is smoother as compared to the Newtonian 
one, probably because in ND momentum conservation 
law require many body interactions for cluster aggrega- 
tion. Two body interactions between two monomers can 
not produce a bounded dimer state. Once a small number 
of small clusters are present in the system, aggregation 
speeds up significantly. 

Fig. |2t reports, for ND, the time dependence of the en- 
ergy for the quench at Tf — 0.15. Comparing Fig. |2K 
and Fig. we notice that the major difference with the 
Tf = 0.05 case is in the long time scale. In the present 
case, instead of remaining frozen, the energy continues 
to drift for all densities we considered (between <p = 0.05 
and <p = 0.50) indicating that the separation process does 
not arrest during the simulated time. 

The ^-dependence of the final values of the energy is 
shown in Fig.|3t for the three cases discussed above. For 
the Tf = 0.05 quench, on the overall range of densities 
this quantity varies between 3.1 and 2.8 and it presents a 
maximum close to the critical packing fraction. In equi- 
librium conditions, one would expect the energy to grow 
with the packing fraction but, in this out-of-equilibrium 
situation, this is note the case. To understand how the 
maximum arises, we show in Fig. [3Jd the distribution 
of bonds for three representative packing fraction, 0.10, 
0.25, 0.50. While the high and the low density cases have 
roughly the same average number of bonds, at the crit- 
ical packing fraction such an average is lower, as indeed 
expected by the value of the energy. This difference is 
originated by the presence of a large tail of particles with 
low number of bonds. In other words, at the critical pack- 
ing fraction the aggregate, possess more surface particles 
than at lower and higher packing. This phenomenon may 
results from the fact that at the critical packing fraction 
critical fluctuations are stronger then at lower and higher 
packing. 

It is interesting to note that, within numerical error, ND 
and BD results for Tf = 0.05 arc in agreement. Indeed 
this is not a trivial issue, since one would expect micro- 
scopic dynamics to play an important role in an out of 
equilibrium situation. 
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For Tf = 0.15, the final U values are smaller than 
the ones at Tf — 0.05, providing further evidence that 
at low T a dynamic arrest transition has taken place, 
blocking the approach to an equilibrium phase separated 
state. Indeed, if equilibrium would have been reached, 
U(T = 0.05) should have been lower than U(T = 0.15). 



C. Structure of the aggregates in real space 

An insight of the evolution of the structure during the 
phase separation process can be gained by examining 
Fig. where snapshots of the <fi = 0.10 system are pre- 
sented for different times, for Tf — 0.05 (both for ND 
and BD) as well as for Tf = 0.15. Before structural rear- 
rangements start to take place (top row), the particles are 
homogeneously distributed within the simulation box. In 
the second set of snapshot, the phase separation pro- 
cess starts to create regions richer and lower in colloidal 
particle concentration (second row), which progressively 
coarsen (third row). In the final configurations (forth 
row), both the ND and BD results present a very sim- 
ilar ramified structure for the Tf = 0.05 quench whilst, 
as expected, the structure for Tf = 0.15 is more com- 
pact. Comparing these intermediate states with the last 
recorded one, it appears that for Tf — 0.05 the system 
dynamics has arrested. For Tf — 0.15, however, even 
if the system forms a percolating space-spanning struc- 
ture, its evolution is not arrested. Indeed the system is 
still able to rearrange itself, suggesting the presence of a 
significant mobility. 

It is important to stress that the arrested structures 
obtained by ND and BD are indeed very similar, 
which provides a support to the independence of the 
non-equilibrium structure to the microscopic dynamics. 
Thus, if one is interested, as in the present work, to char- 
acterize the arrested structure, it is possible to use equiv- 
alently BD or ND since the results is insensitive to the 
evolution scheme chosen. Obviously, since BD evolves 
toward the final structure roughly two order of magni- 
tude slower then ND, the use of ND offers a net gain in 
terms of computational time. If one, however, is inter- 
ested in the kinetic of the process, i.e. how the gel struc- 
ture forms, issues concerning the microscopic dynamics 
should be properly addressed. 

Fig- EJi and Fig.^J) show snapshots of the final configura- 
tions for packing fractions between 0.01 and 0.10 for the 
two quenches. At low temperature, i.e. Tf = 0.05, and 
low packing fraction, i.e. (f> — 0.01, the system is made 
up of large elongated cluster that does not span all the 
simulation box, i.e. does not percolate. Percolation is 
not observed at any stage during the separation process. 
Some monomers are still present, but their number pro- 
gressively decreases with time. In real system the large 
aggregates will be eventually subject to precipitation due 
to the gravitational field. At cf> — 0.05, the final structure 
percolates. All particles belong to the spanning cluster, 
which is extremely open and with tiny connections. The 



situation is even more evident at <fi = 0.10. The net- 
work that span the box has a more defined structure 
and distribution of particle appears to be more uniform. 
This structure clearly resemble the idea of a gel, highly 
inhomogeneous percolating arrested structure. It is per- 
haps interesting to stress that the percolating structure 
is formed during the separation process. 
For Tf = 0.15, at low packing fraction, i.e. (f> = 0.01, 
the final structure is similar to the previous case, the 
system has formed a single non percolating cluster. Its 
structure, however, is more compact than in the previ- 
ous case resembling more the shape of spherical droplet. 
At 4> — 0.05, the situation is different. A single cluster 
is present but it does not form anymore a rigid perco- 
lating structure but rather an elongated cluster. It is 
only aX 4>> 0.10, that the systems percolates again. The 
general trend however, is that the structures are signif- 
icantly more compact than for Tf = 0.05, in agreement 
with the larger number of bonds (lower energy) observed 
at Tf = 0.15 as compared to the Tf — 0.05 case. 

D. Static structure factors 

To properly quantify the structure of the system we 
study the static structure factor as a function of the time 
elapsed after the quench. For a binary mixture the partial 
static structure factors, S a p{q,t) can be defined as 

S afi (q,t) = (fl*(q,i)fi /J (q 1 t)) (4) 

where the partial density variables are defined as g a (q) = 

J2k=i ex P[*Q ' r (*)i Q ^]/\ // A. At equilibrium, time trans- 
lation invariance implies S a p{q,t) = S a p(q,0), i.e. the 
static structure factor is a time independent function. In 
our case, however, we are in an out of equilibrium condi- 
tion and the t dependence must be maintained. In par- 
ticular we shall focus on the total static structure factor 
defined as: 

1.2 

S(q,t)=j2S af3 (q,t) (5) 

Numerically this quantity is calculated for a given con- 
figuration at a time t by calculating the density variable 
Qa (q) and averaging on different q- vectors with the same 
modulus q but different orientation. Such spherical av- 
erage was performed over up to 300 distinct g-vectors. 
Results for the ND Tf — 0.05 case are shown in Fig. 
for three different densities. ((/> = 0.10 in Fig.^J,), critical 
{(f) = 0.25 in Fig.Ep) and high (0 = 0.50 in Fig.Et). 
At <f> — 0.10, the structure factor starts from the typical 
shape of a low density hard sphere fluid, i.e. no structural 
peaks either at contact, i.e. q ~ 6.2 or at higher values. 
As the energy starts to drop, roughly at t ~ 10 2 , the 
system becomes more structured. A first effect is repre- 
sented by an increase at the nearest neighbor peak of the 
structure factor, a structural confirmation of bond for- 
mation. Simultaneously, a spinodal decomposition peak 



5 



(at low q) starts to emerge. This indicates the develop- 
ment of a characteristic length roughly proportional to 
2n/q p , where q p is the location of the peak. This peak 
grows, coarsen and eventually stops. 
To characterize the evolution of the peak we studied the 
evolution of the maximum and of the first moment q\ of 
S(q). This moment is defined as: 

qi = fgc g , , w (6) 
J S(q,t)dq 

where the integral is performed with a proper cutoff q c . 
This quantity scales as the peak location but can be cal- 
culated with better accuracy [42], |5fJ ■ We have chosen 
q c — 4.0, to weight the long wavelength fluctuations re- 
sponsible for the phase separation. S(q p ,t) and qi{i) are 
shown in the inset of Fig.|SJi. The maximum of the peak 
start to grow at t ~ 10 2 and then at t ~ 4 x 10 3 reaches 
a value around 20 and it remains frozen at this value. 
This is also confirmed by the behaviour of q\ that af- 
ter a transient time reaches a plateau. As discuss before 
[39l Hot Hg | , these are phenomena that are typical of the 
formation of a dynamic arrested structure [3l|, • 
For density close to the critical one, i.e. <fi — 0.25, the 
behaviour is similar to the previous case, but less pro- 
nounced, as shown in Fig. HJJd. S(q,0) presents a more 
pronounced contact peak, a consequence of the higher 
density. As the system evolves in time the contact peak 
grows as in the cj> = 0.10 case. The dynamical evolution 
of the maximum and of the first moment qi, shown in 
the inset of Fig. is similar to the previous case but 
the final plateau values are much lower. Hence in this 
situation the arrested phase presents smaller aggregates 
and is more homogeneous then the case at lower density, 
quantifying the previous observations based on the real 
space pictures (see Sec. lIIlCj) . 

A different situation emerges when the system is at high 
density. The evolution of the structure factor for <fr = 0.50 
is presented in Fig. |SJ:. S(q, 0) presents a very pro- 
nounced contact peak, higher then 2.5 due to the high 
particles packing. After the quench, when the system 
evolves, this first peak decreases and moves to slightly 
higher q, i.e. the opposite trend of the previous case. 
This is not surprising, it is well known that when the 
system at high density an temperature is cooled the ef- 
fect for SRAC systems is the lowering of the first peak 
and the appearance of long tail oscillations. Indeed in 
MCT this trend is responsible for the reentrant glass line 
and for the formation of the attractive glass |18j and it 
has been observed in simulation on SW systems both at 
equilibrium |5l| and in the aging regime poj . It is in- 
teresting to note that, while the contact peak decreases 
with time, the second peak (q ~ 12.5) increases, a clear 
indication of the emergence of these oscillations. Also fol- 
lower q, the situation is different from the two densities 
discussed above. A slightly increase of in the structure 
factor is present with time but now we can not talk of a 
real peak formation. The system freezes but its structure 
remains basically homogeneous and no large aggregates 



are formed, as confirmed by the evolution of the maxi- 
mum and of qi reported in the inset of Fig. [U. These 
results suggest that at this packing fraction the system 
is sufficiently homogeneous that the MCT scenario is re- 
covered. The frozen structure is similar to the one of 
an homogeneous liquid and the contact interactions are 
responsible for the structural arrest, enforcing the idea 
that SRAC gels are spatially inhomogeneous attractive 
glasses (in the MCT fashion), where the inhomogeneity 
is build by the phase separation process. 
Analysis of the BD configurations confirms the above sce- 
nario, in agreement with the previous discussion based on 
the potential energy. Fig. 0i, Fig. 03 and Fig. 0: report 
BD results for <j> = 0.10, <t> = 0.25 and = 0.50 re- 
spectively. Results are very similar to that we previously 
discussed for ND. For (f> = 0.10 and = 0.25 a low q 
peak emerges and get arrested at a plateau. The differ- 
ence with ND results is represented by the time it takes 
to get to the maximum value of the low q peak. In the 
BD, the transient is slower as was found in the energy 
evolution. Indeed this phenomenon is more evident for 
4> = 0.10, where in the explored simulation window, only 
the approach to the plateau can be studied. At <fi = 0.50, 
no low-g peak emerges in agreement with the ND results. 
The situation is different for the higher temperature 
quench case, i.e. T f = 0.15. Fig. Hi, Fig. 03 and Fig. 03; 
report S(q,t) for <f> = 0.10,0.25,0.50 respectively. As 
for the Tf — 0.05 case, following the quench the system 
starts to develop inhomogeneities, expressed by the for- 
mation of a growing low q peak in S(q,t). In this case, 
however, no signs of structural arrest are observed. For 
example for = 0.10, Fig. [SJi, S(q,t) develops a max- 
imum at t ~ 10 3 that shows a slow but continuous in- 
crease. Similarly the position of the peak, represented 
by qx, moves to lower q and continues to drift for all the 
time of the simulation, indicating that the aggregates 
are growing in size. At later stage, f~4x 10 4 , the slow 
growth of the maximum presents a sharp increase giving 
evidence that a different regime has been reached. Sim- 
ilar behaviour is observed for (f> = 0.20, Fig. 03, but in 
this case the increase in the maximum and position of the 
first peak is faster as expected since the quench is deeper 
in the phase separating region. Finally for = 0.50, 
Fig. 0:, there is an evident formation of a low q peak 
that, differently from the Tf = 0.05, slowly grows. 
To conclude our analysis, we present S(q) of final config- 
urations at all the densities that have been investigated 
with ND. Fig. 01 shows the Tf = 0.05 case. The first 
thing to notice is that the shape of the structure factors 
does not significantly depend on cj> when q > 5.0, sug- 
gesting a similar local structure. At lower q, however, 
differences are evident. For the lowest density, </> = 0.01 
and <j> = 0.05, the system proceeds with the phase separa- 
tion but for higher density this phenomenon is arrested. 
Moreover, increasing <fi the degree of dis-homogeneity in 
the system decreases continuously up to = 0.50 where 
a low q peak is not detected and the structure is very sim- 
ilar to that of homogeneous liquid. From this analysis it 
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is evident that we can distinguish between two different <p 
regions, (i) At low packing fraction, i.e. 4> < 0-10 the sys- 
tem does not present any arrest in the phase separating 
process, (ii) In the intermediate regime, i.e 0.10 < (f>, the 
system gets arrested in a structure that presents a lower 
degree of dis-homogeneity as the density is increased. 
The cross-over between these two regimes is character- 
ized by the connectivity properties of the resulting struc- 
ture. The onset of a percolating attractive glass structure 
appears to be the condition for global arrest. At small 
packing fractions, glass clusters can still freely diffuse and 
aggregate progressively. Formation of a gel built with a 
diffusion limited cluster mechanism in which the glassy 
clusters created in the phase separation process are the 
renormalized monomers is a possibility which can not be 
excluded, but it would requires the study of a much big- 
ger system, which is out of our present capabilities. This 
issue remains open for further investigations. 
In the Tf = 0.15 case, the low q regime grows indefinitely 
since at this value of the temperature the phase separa- 
tion is taking over and the effect of the low temperature 
affects only the kinetics of this process. 



E. Structure of the aggregates in two dimensional 
slabs 



To better characterized this nearly two dimensional con- 
figuration, we calculate the two dimensional structure 
factor S2d(q) defined as: 

1.2 

a, (3 



with, 



N(z) 



£?a(q) 



^ exp i (q x 



» 



x k + lyVk 



(a) 



(8) 



where the sum is restricted to coordinates belonging to 
the particle that lies within z + A/2 and z — A/2, N(z) 
is the number of such particles, n a (z) the relative con- 
centration of the two species in the slabs. The vector q, 
lies in the plane of the slabs and has the chosen modulus 
whereas the orientation is chosen randomly. For the total 
two dimensional structure factor a definition analogous to 
the one given by Eq.[S] holds. Figll2b reports the results 
for Tf = 0.05 for the three z-values considered together 
with the spherical averaged S(q) considered before. In 
all cases in which the slab dimensions are representative 
of the sample, the three dimensional S(q) coincides with 
the two-dimension one. 



In this section we focus our attention on the proper- 
ties of two dimensional slabs of the simulation box along 
one of the axis. This is something that, in general, is 
achieved in confocal microscopy experiments. Indeed, 
the typical size of a colloidal particle allows for this very 
powerful investigation in real space. Among the impor- 
tant phenomena that has been studied by this technique 
we can mention structural relaxation in hard sphere mix- 
tures |52| , nucleation phenomena |53f | , fluid- fluid interfa- 
cial propertie s 154 a nd indeed aggregation and gelation 
phenomenon |55l l56l 15*^ . 

We divide our simulation box in slabs along the z axis, 
mentioning that in our simulation there is no external 
field acting on the system (such for example gravity) and 
consequently the choice of the height of the slabs is ar- 
bitrary. We take slabs of thickness A = 4, i.e. four 
times the diameter of the B particles. We focus on the 
final configuration at 4> = 0.10 and Tf = 0.05, shown in 
Fig. EI] Fig. II lb shows snapshots from the top of three 
arbitrary slabs. The structure is extremely ramified and 
it is possible to notice long elongated clusters that form 
the basic unit of the arrested structure. As expected, the 
distributions of particles in the three slabs are very simi- 
lar since the phase separating process was arrested before 
a big compact cluster was formed. The situation is differ- 
ent for the Tf = 0.15 case. The three slabs present very 
different distributions, showing the presence of a big ag- 
gregate that results from the tendency of the system to 
form a single spherical droplet, since no arrest is taking 
place. On the contrary in a different slab, only few par- 
ticles are present. 



F. Non ergodicity parameter 

When a system forms an arrested state only few of 
the possible configuration in phase space are actually ex- 
plored, i.e. the system is non-ergodic. This is normally 
detected by the density correlation functions, (f>^{t,t w ) 
defined by: 

(j>f{t,t w ) = (g*M,t)g^q,t + t w )) (9) 

where configuration at a time t w after the quench is cor- 
related with a configuration at tw + t. In what follows 
we focus our attention on the total correlator defined, in 
a similar fashion to Eq. [5] by 

1,2 

<MM™) = J>f(i,U (10) 

When the system is at equilibrium, time translational in- 
variance holds and the correlators are independent from 
the waiting time, i.e. 4>^{t,t w ) = (j>% (t). When the 
system starts to lose ergodicity, correlators do not re- 
lax anymore to zero. This effect is usually measured in 
terms of the non ergodicity parameter f q , defined as the 
long time limit of the correlator, i.e. f q = 4>q{t — > oo). 
This quantity represents the order parameter for the glass 
transition, since when it is zero the system is in an er- 
godic state, when it is finite the system is in a non-ergodic 
state. In equilibrium, within MCT formalism it is possi- 
ble to directly calculate f q from the static structure fac- 
tor and eventually test the theory with experimental or 
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numerical results. We studied both the density-density 
correlation function and its long time value. We focus on 
the case Tf = 0.05 and we chose to analyze the system 
closer to the critical packing fraction, i.e. <f> = 0.25 
In Fig. 113b . the evolution of the correlation function, ex- 
pressed as <f> q (t—t w ), is presented for waiting times longer 
then the time needed to thermalized the velocity degree 
of freedom and for a representative g-vector, qoQ = 20. 
For short waiting time, the correlation functions relax to 
zero. For larger waiting time, a non ergodic contribution 
arises. Eventually, for large enough t wi the correlators 
remains on a plateau that is very closed to unity. This 
behavior is different from the evolution of the aging dy- 
namics for SRAC systems at higher density [2(j. In that 
case, after a low temperature quench, there was no ev- 
ident sign of a plateau even at the largest waiting time 
considered. In other words, the aging dynamics of the 
attractive glass is different from the present one, a fact 
that enforce the idea that the two phenomena have dif- 
ferent origins. 

The evolution of f q (t w ) with waiting time as obtained 
from the long time limit of the correlators is shown in 
Fig. 113b for four representative waiting times. The wave- 
vector dependence of f q (t w ) develops progressively, start- 
ing from small q. The width of the non-ergodicity pa- 
rameter is related to localization length. As time goes 
on, more and more particles aggregate, resulting in a de- 
crease of the average localization length. The fact that 
already at short time the system posses a finite low q non- 
ergodicity parameter is a clear indication of the presence 
of a spanning structure made of few particles. As more 
particles join the percolating cluster the system becomes 
more non ergodic on shorter and shorter length scale. 
This suggests that the localization length progressively 
decreases, following the same pattern recently observed 
in both chemical |58j and thermoreversible [53 gels. 
We have attempted to compare the non ergodic behavior 
observed in our simulations with MCT predictions, using 
as input the S(q) calculated from the simulations. Unfor- 
tunately, as well known, MCT overestimates dynamical 
arrest. In the case of HS, arrest is predicted to take 
place for </> > 0.516 if the Percus-Yevich S(q) is used and 
for </> > 0.525 with Verlet- Weiss correction or, more pre- 
cisely, for (f> > 0.546 is the "exact" S(q) calculated from 
simulations is chosen [6(j. In both cases, these critical 
<j) values are smaller than the experimentally and numer- 
ically detected value of (/) = 0.58. Similarly, the MCT 
predictions for the attractive glass overestimate the ideal 
glass transition temperature by more than a factor of two. 
Comparing experimental or simulation data with MCT 
predictions for the ideal glass transition locus requires 
an appropriate mapping in the <p-T plane 24, 61]. MCT 
predictions which do not account for the mapping sug- 
gest that the ideal glass line preempts the spinodal line 
(i.e. it is located above the phase separation curve) [62| . 
Only when the appropriate mapping is accounted for the 
attractive glass line correctly ends in the high colloid con- 
centration side of the spinodal curve 01 ■ It is very un- 



fortunate that the overestimate of the glass critical line 
prevents the possibility of meaningful solving the MCT 
equations for the case of the phase separating system by 
using as input the "exact" numerical S(q,t w ), since the 
role of the large q is dominating already at short time. 
Our attempts to solve the MCT equations failed in re- 
producing the small width observed in f q at short times 
(t = 57 in Fig. 113b). The first non vanishing f q was char- 
acterized by a width already larger than 50 (in units of 
q°B )• 

We perform a further analysis considering the varia- 
tion of the final f q with density. The result is shown 
in Fig. 114b for a packing fraction ranging from 0.05 to 
0.50. The width of the non ergodicity parameter de- 
creases when the packing fraction is decreased. This is 
an indication of the fact that the average localization 
length is larger for the lower packing fraction, the one 
characterized by a more open structure. As we discussed 
above, at high density the system maintains a certain ho- 
mogeneity and no open empty region are detected. An 
open structure presents large amplitude modes and, as 
a consequence, the mean localization length is expected 
to grow with decreasing packing fraction. This inter- 
pretation is confirmed by the mean square displacement 
(MSD). In out of equilibrium, this quantity is defined as 
(\r(t — t w ) — r(t w )\ 2 ), where r(t) is the position of the par- 
ticle at time t and the average (•) is performed over all the 
particles. Here we focus on the total MSD, i.e. evaluated 
with no distinction between particles of the two species 
for different values of t + w. Results are presented in 
Fig. 114b . After a ballistic short time region, the MSD 
reaches a plateau, whose value is larger the smaller <j) is. 
Since the long-time value of the MSD provides an esti- 
mate of the characteristic size of the cages confining the 
particles, data in Fig.. 114b confirm that, in the arrested 
state, particles are more and more localized on increasing 

<f>. 



IV. CONCLUSIONS 

The aim of this manuscript is to provide numeri- 
cal evidence that, in the case of SRAC systems, the 
formation of an arrested state at low packing fraction 
results from a phase separation process interrupted by 
an attractive glass transition. We have shown that 
indeed, only when the system is quenched below T x , the 
temperature at which the attractive glass line crosses 
the coexistence line on the high-polymer concentration 
side, the coarsening dynamics get arrested due to the 
mobility reduction associated to the glass transition. 
We also observed that the arrest takes place only when 
the structure is percolating, which at the studied Tf, 
requires <fi > 0.05. For smaller </>, diffusion of the 
droplets provides a slow coarsening mechanism which 
is missing in the percolating case. For T > T x , we 
observed a continuous progressive coarsening process. 
The percolating ramified structures produced during 
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the early stages of the phase separation process become 
thicker and thicker and the percolating cluster eventually 
collapses in the attempt of minimizing the surface area. 
These slow rearrangement processes are possible due to 
the residual mobility of the high <f> phase. 
The final structure of the aggregates as a function of 
<p is characterized by a peak at a finite wave- vector q pi 
which constitute the frozen memory of the interrupted 
phase separation process. The location of the peak and 
its amplitude depend on <f>. The smaller the <f>, the lower 
is q p and the larger is S(q p ). A complex interaction 
between the strength of the phase separation and the 
mobility in the dense phase define the final structure 
of the system. For example, for <f> = 0.50, the system 
is so dense that it freezes before it can get significantly 
inhomogeneous and no considerable peak is detected 
in the static structure factor. It is also interesting to 
note that in this case, the evolution of S(q) is not very 
different from the case of a quench in the attractive glass 
phase region |2jj, i.e. a decrease in the contact peak 
with waiting time. In the low density case the trend 
is the opposite, since more particles experience contact 
after the quench. It is also worth recalling that the 
connectivity of the frozen structure at particle level is 
also a function of <fi. We find that the number of bonds 
is minimized roughly in correspondence of the critical 
packing fraction. Perhaps this is related to the strength 
of the critical fluctuation, which increase the surface of 
the aggregate lowering, as a consequence the number of 
bonds. 

Finally we studied the evolution of the non ergodicity 
parameter evaluated from the density-density correlation 
functions. We found that, differently from the aging 
dynamics in glasses, the non ergodicity parameter 
progressively increases with t w . The first components 



to become significantly non ergodic are the large wave- 
length density fluctuations. The width of f q increases 
progressively during the coarsening dynamics. When 
dynamic arrest is completed, the localization length 
has become extremely narrow. Moreover it changes 
with packing fraction, showing less localized aggregates 
at lower density. The progressive increase of the f q 
width is reminiscent of the behavior recently observed 
in chemical (5^ and thermoreversible [H^ gel formation. 
A final remark is on the effect of the microscopic 
dynamics. In order to exclude any artifact introduced 
by Newtonian dynamics, we performed simulations, for 
the lowest T, also using Brownian dynamics. The two 
microscopic dynamics provide the same equilibrium 
description but different time scale. On approaching 
a structural glass transition, it has been shown that 
ND and BD generate the same long time behavior |63j . 
When the system is phase separating, as in our case, the 
microscopic dynamics could play a major role. We show, 
however, that the arrested structure obtained by the two 
different dynamical schemes arc very similar. Since BD 
simulation can be a few order of magnitude slower than 
ND simulation, it is important to have the possibility to 
use the latter to describe aging processes in colloids. 
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FIG. 1: Graphical representation of the studied state points 
in the temperature-packing fraction plane. Crosses represent 
the equilibrium starting configurations at T = 1, which, at 
time zero, are quenched at Tf = 0.15 (circles) and Tf — 0.05 
(squares) respectively. Symbols are shaded if, at the end of 
the simulation, the structure of the system is percolating. The 
error bars, calculated as explained in the text, provide an in- 
dication of the spinodal curve for the studied model. The 
present estimates are consistent with the phase coexistence 
calculations of Miller and Frenkel |4S| for the Baxter model 
(dashed line), properly transformed using the virial mapping 
(Eq|3J. The red-star indicates the location of the critical 
point, from Ref. |4fi|. To better frame the location of the 
studied state points, the Figure also shows three calculated 
isodiffusivity curves, i.e. the locus of points at which the dif- 
fusion coefficient D/D = 5 x 10 _2 ,1 x 10" 2 and 5 x 10~ 3 
with D — aB\/~(T/M). The bold line is the extrapolation 
to D/Do — ► from Ref. [42| . rescaled to the present well 
width case using EqE] The extrapolated D / D — > line 
(see Ref.|24j|') provides an estimate of the location of line of 
dynamical arrest for this system. 
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FIG. 2: (a) Time evolution of the potential energy per particle 
U/N at Tf = 0.05, following the quench from T = 1, for 
different packing fractions, (b) Same as (a) but for Brownian 
dynamics, (c) Same as (a) but for Tf = 0.15. 
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FIG. 3: (a) Potential energy per particle at the end of the 
simulation versus <f> for Tf = 0.05 (Newtonian and Brownian 
dynamics) and for Tf = 0.15 (Newtonian), (b) Histograms 
of the distribution of bonds per particle for three packing 
fractions: 6 = 0.10, d> = 0.25, 6 = 0.50. 



FIG. 4: Snapshots of the configurations at different times 
from the quench at Tf — 0.05 and <j> = 0.10: (a) Newtonian 
Dynamics and (b) Brownian Dynamics, (c) same as (a) for 
T f = 0.15 and (j> = 0.10 



(T/ = 0.05) 



(T f = 0.15) 



FIG. 5: Snapshots of the final configurations for Tf = 0.05 
(top left) and Tf = 0.15 (top right) for 0.01 < 4> < 0.5. The 
unity of length (at) is graphically kept constant for the first 
and the last three rows. Hence, an increase in density implies 
a smaller simulation box. 
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FIG. 6: Time evolution of the static structure factor S(q) 
for Tf = 0.05 at three packing fractions: (a) (f> = 0.10, (b) 
4> = 0.25 (the critical packing fraction), (c) <j> = 0.50. In 
the inset evolution of the height of the first maximum and its 
position expressed by the quantity qi is shown(see text for 
details). The full symbols in the inset represent the times at 
which the structure factors are shown in the main figure. 



10 




5 10 ' 15 



FIG. 9: Static structure factors S(q) of the final configuration 
for different packing fractions, (a) Tj — 0.05, (b) Tf = 0.15 



FIG. 10: Representation of the studied two dimensional slabs. 
The three slabs are chosen at an height z = 2.4, z = 12.0 and 
22.6 and they have a width of 4<T{, 



FIG. 11: Snapshots of the two dimensional slabs (see Fig llOU 
of the final configurations at three values of the ^-coordinate, 
(a) T f = 0.05 and cf> = 0.10; (b) T f = 0.15 and = 0.10 




FIG. 12: Two dimensional structure factor S , 2d(g) for dif- 
ferent values of the ^-coordinate (see text for details), (a) 
Tf = 0.05 and <j> = 0.10; (b) T f = 0.15 and <j> = 0.10. The 
dashed line corresponds to the three dimensional spherically 
averaged S(q). 
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FIG. 13: a) Density-density correlation function 4> q (t) for 
different waiting times t w after the quench. This case cor- 
responds to (f> = 0.25 and T/ = 0.05. Correlation functions 
are calculated for qoB = 20. (b) Time evolution of the non 
ergodicity parameter f q with increasing t w for the same case. 
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FIG. 14: a) Non ergodicity parameter for different packing 
fraction in the final arrested state, at Tf = 0.05. 
b) Mean square displacement in the final arrested state for 
different values of d>. 



